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ABSTRACT 

We investigate the structure of the dark matter halo formed in the cold 
dark matter scenarios by N-body simulations with parallel treecode on GRAPE 
cluster systems. We simulated 8 halos with the mass of 4.4 x lO^^M© to 
1.6 X IO^^Mq in the SCDM and LCDM model using up to 30 million particles. 
With the resolution of our simulations, the density profile is reliable down to 
0.2 percent of the virial radius. Our results show that the slope of inner cusp 
within 1 percent virial radius is shallower than —1.5, and the radius where 
the shallowing starts exhibits run-to-run variation, which means the innermost 
profile is not universal. 

Subject headings: cosmology: theory — dark matter — galaxies: clusters: general 
— methods: N-body simulations 

1. Introduction 

Since the "finding" of the universal profile by Navarro, Frenk, and White (1996, 
1997, hereafter NFW), the structure of dark matter halos formed through dissipationless 
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hierarchical clustering from cosmological initial setting has been explored by many 
researchers. NFW performed a number of A^-body simulations of the halo formation using 
10-20k particles and found that the profile of dark matter halo could be fitted to a simple 
formula (hereafter, NFW profile) 

(1) 



(r/ro)(l + r/ro)2 

where po is a characteristic density and Vq is a scale radius. They also argued that the 
profile has the same shape, independent of the halo mass, the power spectrum of the initial 
density fiuctuation or other cosmological parameters. 

Several groups reported the results of similar simulations with much higher resolutions. 

However, disagreement concerning the inner structure still remains. Some researchers 
claimed that the slope of the inner cusp is steeper than that in the NFW's results. Fukushige 
and Makino (1997) performed a similar simulation with 768k particles and found that the 
galaxy-sized halo has a cusp steeper than p oc r^^. Moore et al. (1998, 1999 hereafter M99) 
and Ghigna et al. (2000) performed simulations with up to 4M particles and obtained the 
results that the profile has a cusp proportional to r~^'^ both in galaxy-sized and cluster-sized 
halos. M99 proposed the modified universal profile (hereafter, M99 profile), 

(r/roMllir/roY-^Y 

Fukushige and Makino (2001, Paper I; 2003, Paper II) performed two series of iV-body 
simulations, and found that the halos have density cusps proportional to r~^^^, independent 
of the halo mass and cosmological models. 

On the other hands, other researchers obtained the slope of inner cusp shallower than 
— 1.5 and close to that in the NFW profile. Jing and Suto (2000, 2002) performed a series of 
A^-body simulations and concluded that the power of the cusp depends on mass. It varies 
from —1.5 for galaxy mass halo to —1.1 for cluster mass halo. Klypin et al. (2001) obtained 
the slope at the center that could be approximated by r"^-^. They, however, argued that 
the NFW fit is still good up to their resolution limit. Power et al. (2003) simulated an 
LCDM galaxy-sized halo with 3M particles and claimed that their circular velocity profile 
obtained is in better agreement with the NFW profile than with the M99 profile. 

The purpose of this paper is to explore the inner structure of the dark matter halo 
by means of A^-body simulations with about 10 times higher mass resolution than that of 
previous simulations. We simulated the formation of 8 cluster-sized halos in the SCDM and 
LCDM models using parallel Barnes-Hut treecode on parallel GRAPE cluster. 



The structure of this paper is as follows. In section 2, we describe the model of our 
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A^-body simulation. In section 3, we present the results of simulation. Section 4 is for 
conclusion and discussion. 



2. Simulation Method 

We consider two cosmological models, SCDM model {Qq = 1.0, h = 0.5, = 0.6) and 
LCDM model {Qq = 0.3, Aq = 0.7, h = 0.7, as = 1.0). Here, Qq is the density parameter, 
Ao is the dimensionless cosmological constant, and Hq = 100/ikm/s ■ Mpc~^ at the present 
epoch. The amplitudes of the power spectrum in CDM models are normalized using the 
top-hat filtered mass variance at 8h~^ Mpc according to the cluster abundance (Kitayama 
& Suto 1997). 

We simulate the formation of the dark matter halos using the "re-simulation" method, 
which has been a standard method for the simulation of halo formation since NFW (1996). 
The procedure for setting the initial condition of halos are the same as that used in Paper 
II. 

We first performed large scale cosmological simulations with 3.7 x 10^ particles in 
a sphere of 300/iMpc comoving radius. We regard spherical overdensity regions around 
local potential minima within as candidate halos. We define the radius such that 
the spherical overdensity inside is 178f2Q^ times the critical density for SCDM and 178f^o^ 
times for LCDM model (Eke, Cole, Frenk 1998). 

We selected 8 halos from the catalog of candidate halos. The selected halos are 
summarized in Table 0. We selected the three most massive halos and one halo randomly 
from halo candidates lying within 200/;,Mpc from the center in both models (so that the 
external tidal field can be included). We express a region within 5r^ from the center of 
the halo at 2; = in the cosmological simulation with larger number of particles. We place 
particles whose mass is as same as that in the cosmological simulation in a sphere of ~ 
lOO/iMpc comoving radius surrounding the high resolution region, in order to express the 
external tidal field. The total number of particles, A^, is listed in Table |1|. The generation 
of initial density fluctuation were done on HITACHI SR8000 (1 node) in Information 
Technology Center, University of Tokyo using C0SMICS2 package (Bertschinger 2001). 

We integrate the system directly in the physical coordinates for both the cosmological 
and halo simulations. We used a leap-flog integrator with shared and constant timestep. 
The step size for the cosmological simulation is At/(tH — ^i) = 1/1024 and that for the halo 
simulation is 1/4096. Here, tn is the Hubble time and ti is the time at which the simulation 
starts. The gravitational softening is constant in the physical coordinates and the length 
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egrav is 5kpc for the cosmological simulation, and Ikpc for Runs SI, S2, LI, and L2, and 
2kpc for other runs of halo simulations. 

The force calculation is done with the parallel Barnes-Hut tree code on GRAPE clusters 
(Kawai, Makino 2003)0. GRAPE is a special-purpose computer designed to accelerate 
A^-body simulations. The parallelization scheme we used are basically the same as Warren 
& Salmon's (1993) Hashed Oct-Tree algorithm, except that we incorporated Barnes' (1990) 
modified algorithm. The modification is necessary in order to make GRAPE work efficiently 
(Makino 1991). We use only the dipole expansion and the opening parameter 9 = 0.4 for 
the cosmological simulation and 6 = 0.5 for the halo simulation. 

For high-resolution halo simulations, we used both a parallel GRAPE-5 cluster at 
University of Tokyo and a parallel MDGRAPE-2 cluster at RIKEN. The parallel GRAPE-5 
cluster consists of 8 host computers (Pentium 4/1.9GHz, 1845) each of which has one 
GRAPE-5 (Kawai et al. 2000) board. The parallel MDGRAPE-2 cluster consists of 8 host 
computers (Pentium4/2.2GHz, 1850) each of which has one MDGRAPE-2 (Susukita et al. 
2003) board. For cosmological simulations we used one board GRAPE-5. The simulation 
presented below needs, for example in Run S2, ~ 300 seconds per timestep, and thus one 
run (4096 timesteps) is completed in 350 hours (wallclock time) with the GRAPE-5 cluster. 

3. Results 
3.1. Snapshots 

Figure ^ shows the particle distribution for Run S2 at z = 0.58 and 0. The length of 
the side for each panel is 6.67 Mpc. For these plots, we shifted the origin of coordinates to 
the position of the potential minimum. In Table |l|, we summarized the radius r^, the mass 
Mv, and the number of particles within aX z = 0. 

3.2. Density Profile 

Figures |^ and § show the density profiles of all runs at 2; = for SCDM and LCDM 
models, respectively. The exception is Run L4, for which we plot the density profile at 
z = 0.06 because the merging process occurs just near the center of halos at 2; = 0. The 
position of the center of the halo was determined using the potential minimum and the 



The source code for both serial and paraUel implementations are available upon request. 
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density was averaged over each spherical shell whose width is log]^Q(Ar) = 0.0172. For the 
illustrative purpose, the densities are shifted vertically. 

We plot the densities by the thick (colored magenta in online edition) lines only if the 
criteria for two-body relaxation introduced in Paper I, tj:e\{r)/t > 3, is satisfied, where ty-ei{r) 
is the local two-body relaxation time given by 

0.065f3 
- GVmln(i?_/.)' 

(cf. Spitzer 1987) and -Rmax is a maximum impact parameter. Here we set Rmax to 1 Mpc as 
a system size. We also confirmed that other numerical artifacts due to the time integration 



did not influence the density profile as will be discussed in section |3.3.2| . The potential 
softening does not influence the profile since the reliability limit obtained by the above 
criterion is more than three times larger than the softening length for all runs. 

At r > 0.02 Mpc or r > O.Olrv, the density profiles are in good agreement to the 
profile given by equation (0) (the M99 profile) in all runs. This result is consistent with the 
previous simulations performed with a few million particles (M99, Ghigna et al. 2000, Paper 
I, Paper II). The fitting here was done using and the least square fit of (p — Pm99)/pm99 
at 0.03 < r < 0.5 Mpc. The scale radii tq obtained by the fitting are summarized in Table 

B 

On the other hands, at r < O.OItv, we can see the shallowing of the cusp from the power 
— 1.5 for all run. The degree of the shallowing seems to increase as the radius decreases, 
which seemingly suggests that the inner cusp profile does not converge to a single slope. 
Moreover, the point where the profile starts to depart from the r~^'^ cusp is different for 
different runs. For example, in Run SI the departure starts at ~ 0.005rv, while at ~ 0.02rv 
in Run L3. This means that the density profile is not universal. 

In Figure ^ we plot the density profiles for all runs scaled by tq and po, together with 
the M99 profile. We can see that at r/r^ < 0.05 all profiles are systematically shallower 
than the r~^'^ cusp, and that in this region run-to-run variation of the profile is significant. 
On the other hands, at r/rg > 0.05 the profiles are in good agreement with M99 profile. 
Although there are some dispersion from M99 profile at r/ro > 0.3, they are not systematic. 



3.3. Reliability 
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3.3.1. Two-body relaxation 



We test the reliability of the criterion (|^)using the simulations of the same initial 
condition as used in Run SI but with several different values for the total number of 
particles (A^). Except for N, we used the same simulation parameters as in Run SI. Figure 
P show the cumulative mass Mr(r) within the radii, 0.1, 0.03, 0.01, 0.005, and 0.003 Mpc, 
as a function of time, for three simulations with 29(Run SI), 14 and 1 million particles 
within Tv. Figure ^ shows the final density profiles for three simulations. The vertical bars 
indicates the reliability limit obtained by the criterion 

In Figure ^ we can see that the cumulative mass evolution obtained in the simulations 
with 29 and 14 millions particle are in good agreement for r > 0.01 Mpc. This agreement 
indicates that our criterion (0.009 Mpc for 14 millions particle run) gives a good reliability 
limit. In Figure ^ we can also see that the density p obtained in the simulations with 29 
and 14 millions particle are in good agreement outside the reliability limit of 14 millions 
particle run (0.009 Mpc). The agreement of the averaged density is somewhat worse than 
that of the density. This is because the averaged density is integrated quantity. Any error 
of the density inside the sphere of radius r affects the average density at radius r. 

Recently, Power et al. (2003) proposed another reliability criterion for the two-body 
relaxation, given by 



Although their function form is different from ours and ignores the dependence on potential 
softening (see. Fig 3 of Paper I), it gives reliability limits similar to ours. For example, 
in Run SI, the reliability limit given by their criterion is 0.007, 0.009 and 0.025 Mpc for 
simulations with 29, 14 and 1 millions particles. These values are within 15% of our limit 



If the stepsize for the time integration is too large, it also influences the profile. We 
check whether the stepsize of time integration used in our simulations is small enough 
by performing the simulations from the same initial model as Run SI but with several 
different stepsize (At). Except for At, we used the same simulation parameters as in 
Run SI. Figure | show the cumulative mass within the radii, 0.1, 0.03, 0.01, 0.005, and 
0.003 Mpc, as a function of time for three simulations with At/(tH — ^i) =l/4096(Run 
SI), 1/2048 and 1/1024. Figure show the profile of the density p for three simulations 




shown in Figure |^. 



3.3.2. Time integration 
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at t/(tH ~ ^i) = 0.78125. We plot the profile at this time since, in the simulation with 
^t/itn — ti) =1/2048, merging process occurs near the center of halos at around z = 0. 

In these figures we can see that larger stepsize makes the central profile shallower. The 
density profile outside of 0.007 Mpc converges even by adapting 1/2048. Therefore, we can 
conclude that the stepsize of At/(tH — ^i) =1/4096 did not introduce any numerical artifact. 

Power et al. (2003) investigated influences of the large stepsize on the profile, and 
showed that the influence depends also on the softening length. They found that if potential 
softening length is larger than an optimal length, e ~ Ar^/ a reliability limit is given 
by 

Ur) . /At\5/6 



15 



f^V, (5) 



and if softening length is smaller than the optical length more timesteps are required than 
that given by criterion d). 

However, an application of Power et al. (2003) 's criterion to our simulation results 
seems to give unphysically reliability limits. For example, in Run SI, the reliability limit 
given by criterion (^ is 0.016, 0.038 and 0.083 Mpc for simulations with 1/4096, 1/2048 and 
1/1024, respectively. From Figure ^ it is clear that these values are far too large. Although 
we do not fully understand the origin of the difference, such difference is possible. The error 
in time integration is very complicated because it depends not only on the stepsize and 
softening length, as Power et al. (2003) showed, but also on the integration scheme selected 
(ex. variable timestep or not, in comoving space or in physical space). 



3.4. Fitting by NFW profile 

In Figure ITO, we fit the density profiles for all runs to the NFW profile. The fitting here 



was done using and the least square fit of (p — Pnfw)/Pnfw at r < 0.5 Mpc (down to 
the reliability limit). The scale radii tq obtained by the fitting are summarized in Table ^ 
We can see that the NFW profile is not in good agreement with simulation results, except 



for Run L3. Figure show the residual, (p — Pnfw)/Pnfw; together with that for the M99 
profile. The agreement with the NFW profile is not good in all radii, while that with the 
M99 profile is not good only in inner region (r < 0.03Mpc). Moreover, we can see that the 
sign of the residuals for NFW profile systematically change, which means the central slope 
of the NFW profile is too shallower. 
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3.5. Evolution 



Figures |T2| and |T3| sliow tlie growtli of tlie density profile for all runs. The virial radii 
and the masses within the virial radius at the redshift plotted are summarized in Table |^. 
We fit these profiles to the M99 profile. The fitting procedure is as same as that for Figure 
^. The scale radii tq obtained by the fitting are summarized in Table |[ 

At the inner region (r < 0.03 Mpc), we can see the density keeps almost unchanged 
from relatively higher redshift for all runs. This fact also can be seen in the evolution of 
the cumulative mass shown in Figure ^. This means that the density at the inner region is 
determined by that of the smaller halo that collapsed at higher redshift. 

The density profile of the outer region is formed as the halo grows and shows 
universality. Moreover, the agreement with the M99 profile at higher redshift is very well 
down to the radius at which the cusp shallowing can be seen at 2; = 0, independent of the 



cosmological model we simulated in this paper. Figure shows the relation between the 
scale radius tq and density po obtained by the fitting. We can see clearly an evolutionary 
pass along a line, po oc r^-^'^, also independent of the cosmological model. 



3.6. Different Fitting 

In section we see that the agreement with the M99 profile is not good at the inner 
region (r < 0.02Mpc), and also that with the NFW profile is worse in the whole range of 
profile in section Therefore, it is worthwhile to fit the results to other profiles. Here, 
we try to fit the results to two different profiles. 

Firstly, we fit the results to a profile that has an inner cusp shallower than that of the 
M99 profile and steeper than the NFW profile [fitting (1)], given as 

^""^ {r/roY [1 + (r/ro)3-] ^ > 



where a is the power of the inner cusp and we set a = 1.3. In Figure |15|, we fit the density 
profiles to the profile given by equation (j^). The fitting here was done using and the 
least square fit of (p — Pni)/Pni for r < 0.5 Mpc (down to the rehability limit). The scale 
radii tq obtained by the fitting are summarized in Table |^. Figure |16| shows the residual 
(p — Pni)/Pni- The agreement is better than both the M99 and NFW profiles. 

We also tried to add another power law region (oc r^) to the M99 profile [fitting (2)], 
given as 



Po 

Co{r/r,y [1 + ir/rX-''" [1 + {r/roY-'] 



P°2 - TT^^TTT^^^^TT^^^T^^^^TTT ( 
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where 

l/C^o = (ro/rcf [l + ro/r,]^■^-^ (8) 

Tc is another scale radius. Ahhough this profile includes more parameters to fit, it is based 
on the observation that two different mechanisms might be working in the growth of the 



halo as suggested by the analyses in section 3.5 



In Figure |T^, we fit the density profiles to the profile given by equation (|^. Here, for 
simplicity, we set /3 = for all runs and, therefore, the equation (|^ becomes 

Pn2 = S (9) 

Co[l + (r/r,)]^-^[l + (r/ro)i-5] 



where 



1/Co = [1 + (ro/rc)]'-' (10) 



The fitting here was done using tq obtained in the fitting to the M99 profile and the least 
square fit of (p — Pn2)/Pn2 at r < 0.5 Mpc. The scale radii Vc obtained are summarized in 
Table ||. Figure |l^ shows the residual, (p — Pn2)/Pn2- As a matter of course, agreement is 
better than that for any other profile, since we increased the number of fitting parameters. 

Unfortunately, in the present simulations, the region that we can use to determine 
which fitting formula is more appropriate is not so large. Further studies with simulations 
with higher resolution and larger number of samples would be necessary. 



4. Conclusion and Discussion 

We performed A^-body simulations of dark matter halo formation in SCDM and LCDM 
models. We simulated 8 halos whose mass range is 4.4 x IO^^Mq to 1.6 X W^Mq using up 
to 30 millions particles. 

Our main conclusions are: 

(1) We found that, in all runs, the slope of inner cusp within O.Olr^ is shallower than 
— 1.5, and the radius where the shallowing starts exhibits run-to-run variation, which 
means the profile is not universal. 

(2) We found that the profile is in agreement with the M99 profile for r > O.Olrv, and 
are not in agreement with the NFW profile. We present different fitting formulae to 
describe the whole range of the simulation results. 
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Although we found interesting features in the inner structure of dark matter halo by 
new simulations with much higher resolution, we could not achieve the final understanding 
of the structure. One question remained is whether the CDM halo has a fiat core or not. 
Another question is whether the same shallowing can be seen in the halo of galaxy or dwarf 
galaxy size. The origin of the inner structure is also still unclear. In order to answer these 
questions, we are now planning to perform larger simulations using new GRAPE cluster 
system. 
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Table 1: Run Properties 



Model 


Run 




Tv (Mpc) 


iVv(xl06) 




m (IO^Mq) 


l + Zi 


SCDM 


SI 


1.58 X 10^5 


3.08 


29.2 


60.3 


5.39 


44.2 




S2 


1.21 X 10^5 


2.84 


31.2 


60.7 


3.86 


45.5 




S3 


1.21 X 10^^ 


2.84 


4.5 


10.0 


26.5 


37.9 




S4 


4.47 X 10^^ 


2.03 


6.9 


13.9 


6.46 


43.4 


LCDM 


LI 


9.61 X 10^4 


2.43 


25.2 


62.8 


3.80 


51.1 




L2 


6.96 X 10^^ 


2.20 


26.0 


59.9 


2.67 


52.7 




L3 


6.49 X 10^^ 


2.15 


7.2 


16.7 


9.01 


47.5 




L4 


4.45 X 10^^ 


1.88 


7.8 


13.5 


5.67 


49.4 



Table 2: Fitting Parameters 



Run 


r-o (Mpc) 
for Pm99 


ro (Mpc) 
for Pnfw 


ro (Mpc) 
for 


rc (Mpc) 

for 


rv (Mpc) 


SI 


1.36 


0.41 


0.70 


0.0014 


3.05 


S2 


1.31 


0.40 


0.68 


0.0014 


2.82 


S3 


0.88 


0.42 


0.60 


0.0036 


2.82 


S4 


0.66 


0.29 


0.44 


0.0023 


1.97 


LI 


0.75 


0.31 


0.50 


0.0023 


2.40 


L2 


0.95 


0.33 


0.52 


0.0015 


2.13 


L3 


0.48 


0.23 


0.34 


0.0027 


2.13 


L4 


0.57 


0.26 


0.38 


0.0024 


1.82 
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Table 3: Parameters at higher redshift 



Run 


z 


Tv (Mpc) 


(Mo) 


ro(Mpc) 


SI 


3.3 


0.16 


1.9 X 10^=^ 


0.085 




1.2 


0.65 


1.6 X 10^^ 


0.25 


S2 


3.3 


0.15 


1.6 X 10^^ 


0.056 




1.2 


0.55 


9.7 X 10^3 


0.20 


S3 


1.2 


0.55 


9.8 X 10^3 


0.29 


S4 


1.2 


0.46 


5.5 X 10^3 


0.20 


LI 


4.5 


0.11 


1.4 X 10^3 


0.11 




1.9 


0.42 


1.2 X 10^^ 


0.20 


L2 


4.5 


0.086 


7.3 X 10^2 


0.05 




1.9 


0.40 


7.0 X 10^3 


0.18 


L3 


1.9 


0.32 


5.4 X 10^3 


0.24 


L4 


1.9 


0.31 


4.7 X 10^3 


0.13 




Fig. 1. — Snapplots from Run S2 at 2; = 0.58. The length of the side is equal to 6.67Mpc. 




Fig. 2. — Same as Figure |1|, but at 2; = 0. 
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Fig. 3. — Density profile of the halos for all runs of the SCDM model at 2; = 0. Only the 
densities plotted in the thick lines (colored magenta in online edition) satisfy criterion (^) 
in section |3]^ at r < r^. The labels indicate the run name. The profiles except for Run SI 
are vertically shifted downward by 1, 2, 3 dex for Runs S2, S3, and S4, respectively. The 
vertical bar above the profiles indicate O.Olrv. The solid curves (colored blue) indicate the 
density profile given by equation (Q) (M99 profile). 
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Fig. 4. — Same as Figure but for the LCDM model. 




Fig. 5. — Density profiles for all runs scaled by vq and po (Table |). The solid (colored blue 
in online edition) curves indicates the density profile given by equation (0) (M99 profile). 
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Fig. 6. — Cumulative mass within radii, 0.1, 0.03, 0.01, 0.005, and 0.003 Mpc, as a function 
of time, for three simulations with 29(Run SI, thick lines), 14(intermediate thick, colored 
blue in online edition) and 1 (thin, colored red) million particles. 
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Fig. 7. — Density profiles for three simulations with 29(Run SI, thick lines), 14(intermediate 
thick, colored blue in online edition) and 1 (thin, colored red) million particles. The vertical 
bars indicate the reliability limits obtained by criterion (^ in section p.2| . 
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Fig. 8. — Same as Figure ||, but for for three simulations with At/{tii — t{) =l/4096(Run 
SI, thick hues) , 1/2048 (intermediate thick, colored blue in online edition) and l/1024(thin, 
colored red). 
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Fig. 9. — Density profiles for three simulations with At/(tH — ^i) =1/4096 (Run SI, thick 
lines), 1/2048 (intermediate thick, colored blue in online edition) and 1/1024 (colored red). 
The vertical bars indicate the reliability limits obtained by criterion (§) in section p.2| . 
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Fig. 10. — Density profiles for all runs. The solid curves (colored blue in online edition) 
indicate the density profile given by equation (NFW profile). 
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Fig. 11. — Residuals (p — Pmqq) /pmqq and (p — Pnfw)/pnfw as a function of radius. 
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Fig. 12. — Evolution of density profile for all runs of the SCDM model. The numbers near 
profiles indicate the redshift. The profiles at 2; = are plotted by the thin lines. Only the 
densities plotted in the thick (colored magenta in online edition) lines satisfy criterion (§) 
in section ai r < r^. The solid curves (colored blue) indicate the density profile given by 
equation (|) (M99 profile). 
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Fig. 14. — The scale density po as a function of the scale radius tq (Table |^) at the redshift 
plotted in Figures and |T3|. The dashed hne indicates po oc Tq^'^. 
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Fig. 15. — Density profiles for all runs at the inner region. The solid curves indicate the 
density profiles given by equations (H) (colored green in online edition) and (|^) (colored 
blue) . 
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Fig. 16. — Residuals (p — Pni)/Pni and (p — Pra)/ Pn2 as a function of radius. 
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